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SECTION 1 
INTRODUCTION 


Over the past 15 years, researchers have established the major role in the 
formation of macrosegregation played by fluid flow phenomena during 
solidification (1-6). Calculations of macrosegregation which take into 
account the effects of fluid motion have achieved good agreement with 
experimental results for specific experimental setups (1-8). The objective of 
the effort reported here is to develop a computer model of solidification, 
including the effects of fluid flow phenomena,* with application to the 
experimental low-g solidification program at Marshal 1 Space Flight Center 
(MSFC). This effort is a continuation of the work performed under contract 
NAS8-33573 in which steady-state and unsteady-state models of horizontal 
axisymmetric bidirectional solidification were developed (9-11). These models 
were limited to calculation of the effects of the interdendritic fluid flow on 
the final macrosegregation for an input temperature field and under the 
assumption of no fluid flow in the bulk melt. This report describes progress 
toward removing these restrictions by coupling the effects of thermally 
induced convection in the bulk melt and of heat flow throughout the ingot to 
the macrosegreat ion calculation in the solid/liquid (S/L) region. Thus, the 
rate of development of the solid, a sensitive process parameter, is calculated 
by the model automat ical lyj the only thermal inputs are the heat transfer 
coefficient at the chill and the temperature of the chill. The user does not 
need to measure internal temperatures for each case, an inherently difficult 
task for the small ingots in use at MSFC. The primary effect of the flow in 
the bulk melt is to change the temperature field and, thus, to alter the shape 
of the S/L zone and the local solidification rates. Only recently has a 
calculation been performed with the coupled effect of bulk liquid flow in a 
vertically solidified axisymmetric ingot (8). To date no published model 
includes the effect of maintaining the energy balance throughout the ingot 
coupled to a complete macrosegregation calculation. 

Section 3 describes the coupling of the energy and macrosegregation equations 
in the S/L region, including results from one-dimensional and two-dimensional 
models. Due to subtleties in the nonlinear form of the macrosegregation 
equations required to adequately model the MSFC experiments, the anticipated 
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method of solution was found to be inapplicable, and in fact no existing 
numerical technique could be used to solve the problem. Consequently a new 
numerical technique was developed which is applicable to a broad class of 
moving boundary problems (12). This technique is implemented in a two- 
dimensional model of solidification of a binary alloy in a rectangular mold. 
The model is limited to cases which have little reverse flow in the mushy 

zone, i.e., cases which are far from freckle formation. 

The model is implemented as an interactive program allowing on-line user 

selection of process parameters and alloy to be solidified, as well as the 
dynamic selection of graphical and tabular output. The user interface is 

described in Section 5; it requires minimal training and no user knowledge of 
programming. The code was written to provide flexibility for future 
development, ease of understanding to experienced programmers through 
readability, and efficiency of machine usage. The use of modern numerical 
techniques provides computational efficiency and gives fast interactive 
response. Documentation of the software is in Section 6. 

Because of the effort expended in developing and implementing the new 
numerical technique, there was insufficient time to complete the calculation 
of flow in the bulk liquid. The preliminary analysis for this task is 
presented in Section 4. 
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SECTION 2 
NOMENCLATURE 

C^ composition of the interdendr 1 1 Ic liquid (wt. pet.) 

C $ final local average composition (wt. pet.) 

G PL ,G PSL ,G PSE ,G PE specific heats of the bulk liquid, interdendr i t ic liquid, 

primary phase solid, and eutectic, respectively (J/gm°C) 

heat fluxes, see Eq. (3.2.1) 

mass fluxes, see Eq. (3 .2.4) 

solute mass fluxes, see Eq. (3.2.10) 

weight fraction of the intermetal 1 ic phase 
eutectic and the terminal solid solutions a 
respectively. 

9l»9s>9e volume fraction liquid, primary phase solid, and 

eutectic, respectively 

2 

g gravity vector (cm/s ) 

AHj.. heat of solidification for component i (cal/mol) 

H l ,H<.,H e specific enthalpy of the liquid, primary phase solid, and 

eutectic, respectively (J/gm) 

molar enthalpies of pure component i as solid and liquid, 
respectively (cal/mol) 

H| enthalpy of the intermetal l ic phase (cal/g) 

AH 2gg heat of formation at 298 k (cal/mol) 


H?_ , HT. 
iS’ iL 


in the 
and 3, 


F»G 


F M’ G M 


F S ,G S 


f | ,f a’ f 8 
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J 


k 

k l 

K 

M 


M, 


n l ,n m ,n s 


N 


Y 


'EE 



P 


P 0 

q 

t 




EE 


V = ^X’V 


Jacobian determinant of the coordinate transformation, 
see Eq. (3*2.1) 

equilibrium partition ratio 

thermal conductivity (J/cm*s*°C) 

2 

permeability (cm ) 

number of equations to solve, see Eq. (3.2.31) 
atomic weight of component i (gm/mol) 

horizontal computational grid size in the liquid, mushy, 
and solid regions, respectively 

vertical computational grid size 

unit normals to the eutectic and liquidus interfaces, 
respectively 

2 

pressure (dynes/cm ) 

2 

ambient pressure, p(X,Y), (dynes/cm ) 

2 

conductive heat flux, (J/cm *s) 

time since the beginning of solidification (s) 

temperature, eutectic temperature, temperature at the 
liquidus isotherm, respectively (°C or °K) 

speed of the eutectic isotherm normal to itself (cm/s) 

velocity of the liquid (cm/s) 
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perpendicular distance from the chill (cm) 


dimension of the casting in the x-direction, see Fig. 3.1 
(cm) 

mole fraction of component i 

distance from the "bottom" of the mold, see Fig. 3.1 (cm) 

dimension of the casting in the y-direction (cm) 

2 

thermal diffusivity (cm /s) 

constant for the solid- and liquid-phases which is 
related to the heat of mixing (cal/mol) 

volume expansion coefficient ("^C" 1 ) 

transformed y-coordinate 

viscosity (gm/cm*s) 

2 

kinematic viscosity (cm /s) 

transformed x-coordinate in Section 3> vorticity (S ^) in 
Section M 

2 

stream function (cm /s) 

density of the liquid, primary phase solid, liquid 

■3 

eutectic and solid eutectic, respectively (gm/cm ) 
transformed time (s) 

heat flux vector which includes a conductive flux and a 
convective flux as defined by Eq. (3-l«9)» (J/cm *s) 



SECTION 3 

COUPLED ENERGY CALCULATION - FORMULATION AND ANALYSIS 


The central effort In the work performed under Exhibit "C n is the development 
of a model with the macrosegregation calculation coupled to the energy 
calculation. All previous macrosegregation models rely on input descriptions 
of the thermal field, an undesirable situation both because the measurement of 
the temperature field in a small ingot is an inherently difficult task and 
because the model predictions of macrosegregation are sensitive to temporal 
and spatial variations in the temperature field. 

The goal of a coupled energy calculation proved to be much more difficult to 
accomplish than expected, primarily due to the lack of an existing numerical 
technique capable of handling the nonlinear complexities of the 
macrosegregation problem. However, the significant progress made toward this 
goal was achieved by developing a new numerical technique for solving problems 
of this type (12). The model with a coupled energy-macrosegregation 
calculation obtained by using the new numerical technique is reported in this 
section. 

3.1 Physical Description 

3.1.1 Problem Definition 

Consider the solidification of a binary alloy in an insulated rectangular mold 
with a chill at one end (Fig. 3.1). in the experiments modeled here, 
solidification occurs over a range of temperatures so that three regions 
develop. The region ahead of the liquidus isotherm is entirely liquid. 
Between the liquidus and eutectic isotherms is a region called the mushy zone 
or the S/L zone which comprises the solid dendrites and solute-enriched 
liquid. At the eutectic isotherm a finite amount of liquid solidifies 
isothermally as the eutectic, and the region behind the eutectic Isotherm is 
entirely solid. 

Flow of the interdendr itic liquid in the mushy zone is due to the shrinkage or 
expansion which accompanies solidification and to gravity acting on the 
density gradient within the interdendr itic liquid. It is this flow of solute- 
rich liquid which results in the macrosegregation in the solidified casting 
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or ingot. The relationship between interdendr i t ic fluid flow and 
macrosegregation was first described by Flemings et al. (3*5) for flow due to 
solidification contraction only, and later for flow driven by gravity as well 
as solidification contraction by Hehrabian et al. (1,6) and Kou et al. (2). 



Figure 3-1* Alloy Solidification in an Axisymmetric Rectangular Mold 


The solidification model described below is used to calculate the 
macrosegregation within the solidified casting by solving the equations for 
the flow of interdendr itic liquid in the S/L zone coupled to the energy 
equation describing heat flow throughout the casting. With the exception of 
the addition of the energy equation the mathematical model of the flow in the 
S/L zone is based on the same assumptions and analysis as used in the models 
which were reported previously (9-11)* In the S/L region of the alloy 
solidifying in the mold shown in Fig. 3*1* the flow of the interdendritic 
liquid is two-dimensional because there is no body force normal to the x-y 
plane. It is assumed that there is no movement of the solid; the flow of the 
interdendritic liquid is modeled as flow through a porous medium where the 
volume fraction available for flow is the local volume fraction of liquid. As 
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before, density is a specified function of temperature which also accounts for 
the variation in the composition of the interdendritic liquid. 

Unlike previous models, the isotherms are not necessarily vertical. The 
vertical temperature gradient is caused primarily by convection of the 
interdendritic liquid and secondarily by the vertical variation in the thermal 
conductivity of the metal due to macrosegregation. 

The phase diagram for the binary alloy and the thermal properties of the alloy 
as it solidifies along with the heat transfer coefficient at the chill are the 
input for the model. The temperature, pressure, velocity field, and 
macrosegregation are calculated by the model. 

3 . 1 .2 The Enthalpy of the Alloy and the Energy Equation 

The external boundaries of the temperature field that is modeled are at the 
edges of the casting; heat flow in the mold is not modeled. There are two 
internal moving boundaries, the liquidus and eutectic isotherms, which are the 
locations of discontinuities in the thermal behavior of the material. 
Consequently there is a different energy equation for each of the three 
regions defined by these isotherms, and a separate relationship must be 
provided to describe the heat transfer across each of these moving boundaries. 
The location and rate of movement of the moving boundaries is not known a 
priori; consequently they must be calculated concurrently with the solution of 
the energy equations. 

The energy equation for each region can be written as 

(3.1.D 

where the functional form of pH, the enthalpy per unit volume, and 4>, the heat 
flux, depend on the region. Here the heat flux includes a term for the 

conductive flux and a term for convective flux. 

In each region the enthalpy is an average of the enthalpies of the phases 
present weighted by volume fraction of each phase. Thus, in the liquid region 



with 


(3.1.3) 


h l = H* + C pL (T-T E ) for T > T e ; 


in the S/L zone 

PH - 9 L P L H L t g s P s H s , (3.1.11) 

with defined by Eq. (3.1.3) and 

H S “ c psL and (3.1.5) 

in the solid region 

PH = 9 e p $e H e + (1-g E )p s H s (3.1.6) 

with H $ = C psE (T-T E ) (3.1.7) 

and H e - H* + C pE (T-T E ) . (3-1.8) 


The reference enthalpies, H* and H* as well as the specific heats C pL , C pSL> 
C PSE and C PE are constants which must be specified for each alloy. Values for 
Al-A.5^ Cu and Sn-15% Pb are given in Appendix A. The remaining variables are 
defined in Section 2. The relationship between temperature and enthalpy as 
defined by Eqs. (3. 1 .2) -(3- 1 .8) is illustrated in Fig. 3.2. 



Figure 3"2. Temperature-Enthalpy Relationship for an Alloy 
with a Eutectic 
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In Fig. 3-2, the curvature in the mushy zone is caused by the dependence of 
the specific heat on the local fraction liquid. There is a jump discontinuity 
at the eutectic temperature due to the release of the latent heat of 
solidification of the eutectic phase, and a corner at the liquidus temperature 
due to a discontinuity in the specific heat of the material. 

The form of the heat flux vector depends on the active heat transfer processes 
in each zone. In the liquid region, there is heat flow due to conduction as 
well as convection, 

* - p l H l v - k T (T)VT , (3.1.9) 

where the form of the thermal conductivity k^. is given in Appendix A for each 
alloy, and the remaining variables are defined in Section 2. In the mushy 
zone, the heat flux has a convective term due to the interdendr i t ic fluid flow 
as well as a conductive term, 

* - g L P L H L v - k T (T,g L )VT , (3.1.10) 

where the thermal conductivity takes into account the relative amounts of 
solid and liquid phase material. In the solid region, the heat flux is due to 
conduction only, but the conductivity does depend on the amount of eutectic 
phase present; thus 


t * - k T (T,g E )VT . (3.1.11) 

Boundary conditions must be specified at all external boundaries. The normal 
flux is zero at all insulated walls, 

3T/3* = 0 and 9* *» 0 on x ■ X for 0 i y i Y , and (3.1.12) 

= 0 and “0 on y » 0 and y = Y for 0 < x < X , (3.1*13) 

where Y is the dimension of the casting parallel to the chill face and X is 
the distance from the chill to the center line or to the end wall. At the 
chill there is some interface resistance so that the surface temperature of 


3-5 



the casting varies with time, and the heat flux is parameterized by both the 
heat transfer coefficient and the temperature of the chill: 

* x - - h(T-T ch|)1 ) on x - 0 forO<y<Y . (3.1. H) 

The temperature of the chill, is constant. 

At each moving boundary there are two boundary conditions because the position 
of the boundary must also be determined. At the tiquidus isotherm, the model 
allows no undercooling so that solid formation begins at the liquidus 
temperature and the heat flux through the isotherm is continuous; therefore 

on x » x. (t,y) ! (3*1.15) 

L I 


where is the flux from the mushy side, Is the flux to the liquid side 
and n £L is the local normal to the liquidus isotherm. At the eutectic front 
the difference in heat fluxes across the interface is given by the latent heat 
released upon solidification of the eutectic constituent times the rate of 
movement of the interface normal to itself: 

on x ■ x^(t,y) (3.1.16) 


where (pH)* is the enthalpy at the solid side of the interface and (pH)* Is 
the enthalpy at the S/L side of the interface. 

The boundary conditions at the liquidus are called "implicit" because they do 
not reference the position or rate of movement of the boundary. The 
conditions at the eutectic are called "explicit" because U^, the rate of 
movement of the boundary, appears in the second condition. 


[V ^mI^ZE “ U IE [ ■ W«] 
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3.1.3 I nterdendr i t ic Fluid Flow 

The interdendr itic fluid flow in the S/L region is described by the equations 
for total-mass and solute-mass conservation together with Darcy's Law for flow 
through a porous medium. The analysis leading to the conservation equations 
is based upon the concept of a volume element or cell within the S/L zone 
which is small enough to be treated as a differential element, yet; large 
enough so that the volume fraction solid within it is equal to its local 
average value. This volume element concept is developed in Ref. 3. There is 
no solid phase movement into or out of the cell; there is no diffusive mass 
flux into or out of the cell; and the temperature and liquid density vary only 
differentially across the cell. With these assumptions, conservation of total 
mass within the volume element is written. 

3 

(p $ g$ + p l 9 l ) “ “ v * p l 9 l^ * (3-1.17) 

(See Section 2 for definitions of the symbols.) If solute enters or leaves 
the cell only by liquid convection, if there is no solute diffusion within the 
solid phase, and if the composition of the liquid within the volume element is 
uniform, then conservation of solute In the volume element is written 

"at ^ c s p s 9 s + c l p l 9 l^ “ " 7 c l p l 9 l v * ( 3 . 1 . 18 ) 

where is the local average solid composition defined by 

Ti <es p S «s> ' kc L p s IT • 

Darcy's Law for flow through a porous medium is used to relate the 
interdendr itic fluid velocity directly to pressure. Thus 

v - - jjj- (VP ♦ P L g) . (3.1.20) 

where U is the viscosity of the liquid. The permeability, K, is a local 

function of solidification time and fraction liquid as described in Ref. 11. 
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Mass fluxes at the edges of the S/L zone provide boundary conditions for the 
conservation equations. The top and bottom of the S/L zone are in contact 
with the mold walls so that there is no normal component of flow; then 

Vy ° 0 at y * 0 and y » Y for x £ <_ x 1 x L . (3.1.21) 

As described in Ref. 1, flow at the eutectic isotherm must compensate for 
solidification shrinkage or expansion of the eutectic phase: 

P LE vn EE + ^ P SE ” P LE^ U £E ° ® on x = for 0 < y < Y (3-1.22) 

where n^ is the local unit normal to the eutectic isotherm and is the 

speed of the isotherm in the n^ direction. It is assumed that there is no 
convection in the bulk liquid, so that 

P ° P 0 + P|_9 l (Y-y) at x = x L for 0 < y < Y. (3.1.23) 

In addition, at the liquidus isotherm the volume fraction liquid is given by 

g^ = 1 at x ■ x^ for 0 < y < Y . (3.1. 2^ ) 

During the initial transient, the left boundary of the S/L zone is in contact 
with the chill face of the mold and Eq. (3.1.22) is replaced by the condition 
of zero normal mass flux: 

v^ = 0 at x ■ 0 for 0 £ y £ Y . (3.1.25) 

During the final transient, the right boundary of the S/L zone is at the mold 
centerline or wall and Eq. (3-1.23) is applied. 

In the previous solidification models (Refs. 9-11) a "local solute 
redistribution equation" and a "pressure equation" were derived from the flow 
equations for use in obtaining a numerical solution. In the current model, 
finite differences are applied directly to the equations obtained by 
substituting Eq. (3.1.19) and Eq^ (3.1.20) into Eq. (3.1.17) and Eq. 
(3.1.18). Details of the numerical analysis are presented in Section 3*2. 
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3.1.4 Hacrosegreqat ion 

After solidification is complete at a point (x,y) in the casting, the final 
local average composition is given by 

1 ‘ 9 E 

s J kC L dg $ + P 
g s =° 

The integral term accounts for the solute within the dendritic solid, and the 
second term in the numerator accounts for the solute In the eutectic 
constituent. 

3.2 Numerical Methods 

3.2.1 Introduction 

The description of heat flow during the solidification of an alloy is a moving 
boundary problem with two moving internal boundaries which define a region of 
both solid and liquid phases. The boundaries of this region are associated 
with discontinuities in the functional relationship between temperature and 
enthalpy, and a numerical solution of the energy equation in the presence of 
such discontinuities requires special techniques. This section describes a 
new method of implementing the boundary conditions that allows an accurate and 
smooth determination of the internal boundary locations where this information 
is critical as in predicting the flow field within the S/L region. 

3. 2. 1.1 Problem Description 

A moving boundary problem is an initial value problem described by a system of 
parabolic partial differential equations with one or more boundaries which 
move in a way that depends on the solution of the differential equations (13). 
There are two boundary conditions to be satisfied for each moving boundary 
because its position is unknown a priori. Moving boundary problems are 
nonlinear problems: the position of the boundary must be determined 

simultaneously with the solution of the partial differential equations. 

There are two types of boundary conditions at a moving boundary. The most 
common is the M exp!icit" (13) condition where one or both of the boundary 


SE 9 E C e]/[ p S (1_ 9e ) + p SE 9 e] (3.1.26) 
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conditions explicitly specifies the velocity of the boundary. In a problem 
with an "implicit" condition, neither boundary condition specifies the 
velocity of the boundary. As shown below, the boundary conditions at the 
eutectic isotherm are explicit while those at the Hquidus Isotherm are 
implicit. This poses an especially difficult numerical problem because most 
existing techniques that work for explicit boundary conditions do not apply to 
implicit boundary conditions, and techniques designed for problems with 
implicit boundary conditions rely on analytical manipulations that are not 
feasible for problems as complex as the macrosegregation calculation (14,15). 
The exception to this is the "enthalpy technique" (16,17) which Is discussed 
below. 

There is another aspect of the macrosegregation calculation which makes it a 
particularly difficult moving boundary problem. The solution to the 
macrosegregation equations is especially sensitive to the rate of movement of 
the mushy zone both because the boundary condition at the eutectic is in terms 
of the velocity of that boundary and because the movement of the entire mushy 
zone must be accounted for by the numerical technique used to solve the 
conservation equations. Thus slight inaccuracies, especially fluctuations, in 
the calculated rate of movement of the boundaries cause substantial errors and 
a failure to obtain convergence in the solution of the macrosegregation 
equations. 

3*2.1. 2 Development of a Solution Technique 

In recent years a great variety of numerical techniques have been developed to 
solve moving boundary problems. A review can be found in Ref. 13* Existing 
techniques were bound to be unsuitable because one or both of the following 
features could not be treated adequately: 

a. the position and rate of movement of the boundaries are crucial to 
the calculation and, therefore, must be accurately determined, and 

b. there is no explicit boundary condition at the liquidus interface. 

Techniques for moving boundary problems can be differentiated by their 
treatment of the computational node points. Fixed-grid techniques carry out 
the calculations on a grid that is stationary in the physical domain, using 
special interpolation formulas in the cells containing the moving boundary.. 
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Variable space-grid techniques operate on grids that deform as each subregion 
changes. The Murray-Land is method ( 1 8) and applications of finite element 
methods (19) are examples of variable grid techniques. Coordinate 
transformation techniques actually transform the equations to a domain where 
the moving boundaries are fixed for all time; the calculation is carried out 
on the fixed grid in the transformed domain. The enthalpy method (16,17) 
deserves special attention because it has been successfully applied to solving 
many moving boundary problems. It can be used to solve an integral form of 
the energy equation on a grid that is fixed in the physical domain without 
reference to the locations of the moving boundaries. 

Previous applications of all these techniques except the enthalpy method 
require explicit boundary conditions to follow the boundary movement, and thus 
they are not applicable without modification at the liquidus interface. The 
enthalpy method was tested extensively for use in this problem, but with it, 
the positions and velocities of the moving boundaries could not be determined 
with sufficient accuracy to give sensible results for the macrosegregation 
calculation. Because of these difficulties, it was necessary to devise a new 
method of applying the boundary conditions. 

3.2.2 Discrete Flux Method 

The philosophy of the discrete flux method is to preserve the physical 
conservation relations both analytically and in discrete form. In effect, the 
finite difference equations "simulate" heat flow through the discrete 
computational mesh. The resulting discrete form of the moving boundary 
conditions is explicit whether the original analytical form is explicit or 
imp I icit. 

The approach is first to transform the equations to body-fitted curvilinear 
coordinates, retaining both enthalpy and temperature as dependent variables. 
Then write conserving finite difference approximations for each of the 
interior and boundary coordinates. Finally, obtain conditions which determine 
the movement of the internal boundaries such that the correct discrete heat 
flux relations are maintained at the interfaces. This last step is the one 
that distinguishes the discrete flux method from other techniques. 
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3.2.2. 1 Coordinate Transformation 

Each of the three regions is transformed to a unit square as shown in Fig. 
3.3. The moving boundaries are now included among edges of the squares. A 
computational mesh is set up in each of the three squares with horizontal mesh 
widths A£ s , A5 m , AC l and vertical mesh width An for all regions. The uniform 
vertical mesh spacing insures that points in adjacent regions coincide along 
the edges corresponding to the internal moving boundaries. 



Figure 3-3* The Coordinate Transformation 

Following the methods of Thompson et al. (20) and Viviand (21), the energy 
equation for each region is transformed to the computational space and the 
resulting equation is rewritten in conservation law form. 



where F » -y x pH y ♦ - x $ and (3.2.2) 

7 n t 'n x n y w ' 
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J ■ x^y n *be Jacobian determinant of the transformation and 
Xj., x^, x t and are coefficients of the transformation that depend on the 
size and rate of deformation of the corresponding physical region (16). In 
Eqs. (3.2. 1 )— (3-2-3) » F and G are the heat fluxes in the £- and n- directions, 
respectively. The first term in F is due to the movement of the physical 
domain; x„ is the local rate of movement. 

The mass conservation equations are of the same form as the energy equations 
and they transform in the same way. The total mass conservation equation 


transforms to 


H ■ - 7 • (p lV> ^ £ < Jp > - - ( -5T + 

(3.2.A) 

where p » p $ g s + P L g^ , 

(3.2.5) 

F « - - V * p iA <y n v x ‘ *nV> and 

(3.2.6) 

S N ‘ ‘>l9u (X 5 v Y ) • 

(3.2.7) 


F u and G u are mass fluxes in the 5- and n- directions, respectively, relative 
n n 

to the fixed (C,n) coordinates due both to convection of the interdendr i t ic 
liquid and to deformation of the physical (x,y) domain. The solute mass 
conservation equation, Eq. (3.1.18), is slightly complicated by the 
substitution of Eq. (3.1.19). 

3C 

" p at^ + Jt " " 7 ’ (p l 9 l c l^ (3.2.8) 


where 


p => Pg9 s k and PC ■ P^g L C L + P S 9 S k ^L * 


Eq. (3.2.8) transforms to 

[f ( jc l> ♦ K < F c>] 



(3.2.9) 


P 


( 3 . 2 . 10 ) 



where 


F C " ’ y n x r C L * 


(3.2.11) 


F s “ - Yn x x pc + p L 9 L C L (y n v X " x n v Y ) ’ and 


( 3 . 2 . 12 ) 


G S " p l 9 l c l (x c v y ) 


(3.2.13) 


3. 2. 2. 2 Finite Difference Approximations to the Energy Equation 
In order to write conserving finite difference expressions, the energy 
equation is integrated over each cell in the computational mesh, and Green's 
theorem is applied to give 


// 

cel 1 


g^-( JpiD d£dq 


* 


Fdn - Gd£ 


ceil 

boundary 


(3.2.14) 


Thus, the difference equations can be written in terms of fluxes through the 
cell walls. 

3. 2. 2. 2.1 Interior and Wall Ceils 

The (i,j) interior cell and the ( 1 , j) wall cell are illustrated in Fig. 3.4. 
The following discussion applies to all regions. Note that the heat fluxes 
are defined to be positive in the direction of increasing £ or n. Assuming 
that the function values at the node points represent the function throughout 
the cell and applying Eq. (3.2.14) to the interior cell gives 


i — n+1 — n. , ,n+1 . , n+1 /, „ 

7 — (JpH - JpH )..At»A5 * (F. - F JAt» + (G. - G Ja? . (3-2.15) 

At Ij 'in out in out 


F. and F . are obtained by evaluating Eq. (3.2.2) at i-i and i+i 
in out 

respectively, where for example, 


(* x > i+± f j = “k(1/x ? ) (T. +1 j - T. j) /AS, and (3.2.16) 
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Figure 3-^ - Schematic of Heat FJuxes at the Interior and Wall Cells 

( VlH,j = ' l [' ( 'n / V. )(T i*l.j ‘ T i.J )/4E 

* J(l/y n )( T i jtl - T.j., * T i+1J+1 - T i+tJ .,) /2An] <3- 2 - 


When fully expanded, Eq. (3-2.15) is identical to a space-centered finite 
difference approximation. The equation for the wall cell is similar except 

for the difference in cell size; accordingly 

i . 4 i n n +1 n +1 

^:(JpH n - =(F c - F out )An + ( 6 |n - G out )AE/2 . (3.2.18) 

3. 2. 2. 2. 2 Interface Cells 

There remain to be defined the computational expressions for determining the 
movement of the internal interfaces or moving boundaries. These will be 
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derived by applying the method used for the interior cells and the wall cells 
to the cells situated at the moving internal boundaries and then by using the 
moving boundary conditions to resolve the unknown fluxes. 

Figure 3.5 illustrates cells at the internal moving boundaries of the 
computational mesh. The cell on the left is in the solid region and the cell 
on the right is in the mushy zone. When mapped back to physical space, the 
boundary node points from each cell coincide and lie on the eutectic isotherm. 



Figure 3-5. Schematic of Heat Fluxes at Cells in the Solid and 
Mushy Regions along the Eutectic Interface Boundary. 


Applying Eq. (3.2.14) to the cell in the mapped solid region leads to 


- JpH^) N J 4nAC s /2 


(FT - 
in 


out 


n+1 

)An 


+ (G 


in 


, n+1 , 

- G S JAE /2 
out s 


1(3.2.19) 


where the subscript S means the quantity is evaluated in the solid region. In 
this case is known as a result of the first of Eqs. (3.1.16), and the 
temperature-enthalpy relationship, but F* yt ' s considered to be unknown. 
Similarly, the energy conservation equation applied to the cell in the mushy 
region leads to 


At^ JpH M " JpH M^1,j &nA5 M /2 


(F? - 
in 


n+1 
F" JAt, 
out 


.M 


(G? - G M JAC m /2 , (3.2.20) 

in out “ 


n+1 
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where the subscript M means the quantity is evaluated in the mushy region. 
Again, pH^ is known, but F . ^ is unknown. The second moving boundary 

condition at the eutectic isotherm provides another relationship between the 
unknown fluxes; the relationship is 

y x (pH n+1 - ph!] +1 ) * pS * ‘ 'i (3.2.21) 

’ t \ t s M out in . 

The two unknown fluxes are eliminated algebraically from Eqs. (3.2.19)* 

(3.2.20), and (3.2.21), and the resulting equation is used to determine the 

rate of boundary movement by solving for x T . The same procedure is applied at 

the liquidus interface. In this case x T does not appear in the interface 

condition analogous to Eq. (3.2.21), but it does appear in the remaining S- 
M L 

fluxes, F. and F . so that the final condition at the liquidus boundary is 
explicit even though the original form of the boundary condition in real space 
and time is implicit. 

In effect the rate of the internal boundary movement is determined so as to 
satisfy both the energy conservation equations in the cells adjacent to the 
boundary and the condition at the moving boundary itself. The equations are 
well posed for determining the boundary movement so long as temperature and 
enthalpy are not uniform throughout all adjacent cells. 

3. 2. 2. 3 Finite Difference Approximations to the Mass Conservation Equations 
The derivation of the finite difference approximations for the mass 
conservation equations parallels that for the energy equation. The equations 
are integrated over each computational cell, Green's theorem is applied to the 
divergence terms, and then the finite differences are written in terms of mass 
fluxes through the cell walls. For example, the total mass conservation Eq. 
(3.2.4), written as a difference equation in the interior cell shown at the 
left of Fig. 3.4 is 

ji PlH- <F H(ln) -F H(out) )"V ♦ (ShUM-W))^ ' <3 2 - 22) 

In order to apply finite differences to the mass flux terms, the expressions 
containing velocity in Eq. (3*2.6) and Eq. (3*2.7) must be expanded in terms 
of pressure by substituting Darcy's Law, Eq. (3.1.20); this yields 
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(3.2.23) 


F 


M 


where 


y ri x T p “ P L v { C £S If + C £n In + ^ p L“ p L0^ y n 9 x“ x n 9 y^ } 

C £C " (y n + % )/x C y n and C en " ’ V y n ’ and 

G M = " P L 5 { C 5n ff + C nn “ + (p L“ p L0 )(x 5 g y ) } 


(3.2.24) 

(3.2.25) 


where C = x_/y . 

nn E 7 n 

Then finite differences are substituted for the pressure derivatives: 


(3.2.26) 


iR 

3E 




(3.2.27) 




i , j±i 


1 + p l+1,J.t ■ P l i l,j-l l/24n > 

(3.2.28) 

+ p i+!,j ■ p i-1,j )/2St ’ and 

(3.2.29) 

± (p i,J ± l * p 1,j )/iI ' • 

(3.2.30) 


The remaining terms in Eqs. (3*2.23) and (3.2.25) are evaluated by suitable 
averages. The derivation of the equations for the boundary cells is similar. 
The type of boundary condition applied at each side of the S/L zone depends 
upon the stage of the solidification. During the initial transient, for 
example, the component of velocity normal to the left side is zero, and the 
transformed form of Eq. (3.1.25) must be substituted into the wall-cell form 
of Eq. (3.2.25). After the S/L zone is fully developed, the transformed form 
of Eq. (3.1.22) must be substituted into the wall-cell form of Eq. (3*2.25). 


Finite difference approximations to the solute-mass conservation equation are 
derived in an analogous manner. 
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3.2.3 Solution of the System of Nonlinear Equations 

The result of the analysis in the foregoing sections is a large system of 
nonlinear algebraic equations which give the temperature and enthalpy in all 
three zones as well as the volume fraction of liquid, velocity and pressure in 
the S/L zone, and the shape and rate of movement of the S/L zone. The system 
comprises equations which include the temperature-enthalpy relationships, Eqs. 
(3-1 -2)— (3-1 -8) , the finite difference form of the energy equation in all 
three regions, the moving boundary conditions, Eqs. (3-2.19)~(3.2.21 ) , the 
finite difference form of the total-mass conservation equation and the finite 
difference form of the solute-mass conservation equation. If the horizontal 
dimensions of the finite difference grid in the solid, mushy and liquid 
regions are N^, and N^, respectively, then the number of unknowns to be 

determined at each time step and the number of coupled algebraic equations is 

M = N y j^2(N $ + N m + N l - 1) + 2 N h J (3.2.31) 

where N y is the vertical dimension of the grid. 

The system of equations is solved by the iterative Newton-Raphson method (22) 
for locating the zeros of a system of nonlinear equations. The recursive 
equation at each iteration is 

x (k+1) = x (k) _ j-1 q( X 0<)) (3.2.32) 

where X is the vector of M unknowns, Q is the vector of functions whose zero 
is sought, and J is the Jacobian matrix of Q, with respect to X. 

For even a relatively coarse mesh the system of equations represented by Eq. 
(3.2.32) is quite large. For example, if there is a 10 x 10 grid in each 
region then M *» 780 and an evaluation of Eq. (3.2.32) requires solution of a 
780 x 780 matrix equation at each iteration. However, the structure of the 
equations in Q leads to a sparse partitioned structure for J which can be 
reduced analytically to the solution of a sequence of banded tridiagonal and 
pentadiagonal matrix equations, reducing the number of computer calculations 
immensely. 
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3.3 Program Results 


3.3*1 One-Dimensional Results 

For the purpose of developing and testing the numerical techniques for solving 
the moving boundary problem, a great deal of work was performed with a one- 
dimensional heat flow model. Some of the results are presented here. 

3. 3 . 1.1 Comparison of the Enthalpy Method to the Discrete Flux Method 
Figure 3.6 illustrates the difficulty encountered in applying the enthalpy 
method (16, 17) to a one-dimensional macrosegregation calculation for an 

aluminum-copper alloy. Figures 3-6a and b are the position and rate of 
movement of the liquidus boundary as functions of time, and Fig. 3.6c and d 
are the position and rate of movement of the eutectic boundary as functions of 
time. The dotted lines were calculated by interpolating the temperature 

profile produced by the enthalpy method while the solid lines were calculated 
by the discrete flux method. The cusps in the liquidus time history from the 
enthalpy method are caused by the movement of the "corner" in the enthalpy- 
temperature function over the spatial mesh points. They appeared for all mesh 
spacings tested, up to an order of magnitude smaller than that shown in Fig. 
3.6. The resulting rate of movement of the liquidus fluctuates 
discontinuously about the true value and actually becomes slightly negative at 
about 30 seconds. 

The results for the movement of the eutectic boundary, Fig. 3.6c and d, are 
similar. In this case the fluctuations are caused by the movement of the step 
discontinuity over the nodes in the spatial mesh. Note that if the step in 
the enthalpy-temperature profile is replaced with a sloping line segment by 
the method described in Ref. (23), then there would be two corners near the 
eutectic, each with a behavior similar to that at the liquidus "corner." 

Although the enthalpy method applies spatial finite differences across the 
mesh without regard to the location of the interface, the discontinuities and 
corners in the temperature-enthalpy relationship used to resolve the time- 
difference term of the difference equation are propagated as discontinuities 
in the spatial derivatives of enthalpy and temperature. Thus interpolating 
temperature or enthalpy to determine the location of the interface results in 
the type of temporal fluctuation shown in Fig. 3-6. 
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Figure 3.6. Interface Movement Calculated by the Enthalpy 
Method and by the Discrete Flux Method. 







3. 3. 1.2 Comprison with Empirical Results 

A one-dimensional theoretical result using the discrete flux method is 
compared with the experimental results of Koump et al (24) in Fig. 3*7* The 
heat transfer coefficient and the temperature at the chill were adjusted to 
produce the same average cooling rate at the chill and the same time for the 
liquidus to reach the end of the ingot. In Fig. 3.8 the same comparison is 
made for a slower cooling rate. Here only the temperature of the chill was 
adjusted to match the average cooling rate at the chill. Agreement between the 
calculated results and the measurements is reasonably good considering that 
the model used for the calculations assumes that heat transfer at the chill can 
be characterized by a constant and uniform heat transfer coefficient. In 
addition, perfect agreement between the calculated results and the experimental 
can not be expected because, for the three surfaces not chilled, perfect 
adiabaticity is assumed in the model whereas in the experiments such a boundary 
condition was only approximated. 

3-3.2 Two-Dimensional Results 

Figures 3 -9—3 .16 illustrate calculated results from the two-dimensional model 
used to compute temperature, composition and velocity in A1-4.5S Cu alloy as it 
solidifies. Temperature profiles across the casting at several times during 
the solidification are shown in Fig. 3-9- Any vertical variation in the 
temperature field is imperceptible on the scale of this plot. Figure 3.10 
shows no vertical variation in the position of the eutectic isotherm, but a 
slight vertical slope develops in the liquidus isotherm as it progreses. This 
slope is primarily due to the combined effect of convective heat transfer in 
the S/L region and secondarily due to spatial variations in the thermal 
conductivity of the alloy caused by macrosegregat ion. 

Figure 3-11 shows the final composition of the solid formed by the time the 
run was stopped at t ■ 451 seconds. 

Figures 3-12-3.16 show the movement of the S/L zone across the mold, and 
within the S/L zone the vectors indicate the velocity of the interdendr i t ic 
liquid. The S/L zone is bounded on the left and right by solid lines 
indicating the position of the eutectic and liquidus isotherms, respectively. 
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X(FT) I -n I _ X (FT) 


~ 35°C/MIN FAST COOLING 



TIME (HRS) TIME (HRS) 

EXPERIMENTAL MEASUREMENTS THEORETICAL PREDICTION, h - 0.07 CAL PER 

CM 2 -SEC-°C. TCHILL = 110°C 


igure 3.7. Comparison of the Experimental Measurements of Koump et al (24) 
with Model Prediction for Fast Cooling of A1-4.5$ Cu. 


- 5°C/MIN SLOW COOLING 



PER CM 2 SEC-°C, TCHILL - 528°C 


Figure 3.8. Comparison of the Experimental Measurements of Koump et al (24) 
with Model Prediction for Slow Cooling of AT — 4 . 5% Cu. 
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Movement of the Liquidus and Eutectic Isotherms 
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Figure ,3*13. Final 
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Figure 3*15* Final 




4.50 CU SOLIDIFICATION MODEL !3- JAN-84 08:38:50 


X 



V 


3-31 


Figure 3.16. Final Local Average Composition and interdendr i t ic Fluid Velocity at t - 451 Seconds 


As the completely solid-region grows out from the chill the final local 
average composition is indicated by contour lines at intervals of 0.1 wt. pet. 
Cu. The last figure, in the sequence, Fig. 3.16, is at ^51 seconds from the 
beginning of solidification when the calculation stopped due to failure to 
obtain convergence in the Newton-Raphson iteration procedure. Preliminary 
testing indicated the problem was caused by an inadequate approximation in the 
calculation of the total-mass conservation elements of the Jacobian matrix. 
The maximum reverse flow velocity along the lower boundary is about an order 
of magnitude slower than the isotherm velocity, so freckling is not considered 
to be a contributing factor. An alternative procedure with an uncoupled 
solution of the total -mass conservation equation is partially complete at the 
time of this writing. 
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SECTION 4 

FLOW IN THE BULK LIQUID 


In Section 3 and in the previous reports on this research (9-11), the 

convection in the "bulk liquid" is ignored. The region containing the bulk 

liquid is shown in Fig. 3.1. By ignoring the convection in the bulk liquid, 
the pressure along the liquidus isotherm can be specified as 

p = p 0 + P |_9 (*-y) (4.1) 

where p Q is the pressure at x = X and y = Y (the ambient pressure) and P L is 
the density of the bulk liquid. A more rigorous solution to the 

solidification of an alloy would include a coupling of the convective field in 
the bulk liquid to the convective field in the S/L zone. An approach is 
described in this section. 

A. 1 Physical Description 

In the bulk liquid region shown in Fig. 3-1 there is two-dimensional 
convection because the temperature in this region is not uniform. If the 
thermal properties of the liquid are constant, then the continuity, momentum 
and energy equations (25) are 

V • v = 0, (4.1.1) 

!^ + v*Vv Vp+ W Z v , (4.1.2) 

and 

|I + v • VT - ctV Z T (4.1.3) 

in which v is the kinematic viscosity and a is the thermal diffusivity. 

The assumption of constant properties can be relaxed when solutions to Eqs. 
(4.1.1)— (4.1.3) are obtained numerically. The properties can be assigned 
according to the local temperatures; the approximation neglects nonlinear 
terms which are very small.’ 
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To treat natural convection, Wilkes and Churchill (26) define a pressure which 
the deviation from the initial static pressure; a similar definition (but not 
exactly equal) of pressure is used herein for the pressure in the bulk liquid. 
It is 


P = P - [pg + P Li g (Y-y)] (A. 1.4) 

with P L . as the initial density of the liquid and Pg as the ambient pressure. 
With p defined according to Eq. (4.1.4) it is consistent with the modified 
pressure in Eq. (4.5.1) of Maples and Poirier (9). 


The volume expansion coefficient, defined as 


3p 

e = . 1 k 

P P L 3T » 


is approximated as 


6 °L \ T ■ T , 


(4.1.5) 


Then with g = 0, g » -g and by using Eqs. (4.1.4) and (4.1.5), the 
x y 

components of Eq. (4.1.2) can be written as 


3v 


3v 3v / 3 v 3v 

. x. x 13p^l x x 

3t x 3x y 3y p 3x \ _ 2 _ 2 


(4.1.6) 


3x ay 


and 


3v 


3v 


3v 


+ v — + v + Bg (T-T.) 

3t x ax y ay p ay 3 i 


v ( 


2 2 
3 v 3 v 


2 + 2 
3x 3y 


j (4.1. 


7) 


Ridder (8) considered the natural convection in the bulk liquid above the S/L 
region of an ingot solidifiying upwards. In his analysis, he assumed 3p/3x «* 
0, 3p/3y = 0 and p ^/p ^ = 1 where p is the density of the liquid at a reference 
temperature T. Then 


p L = P L [l-B (T-T)] 
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V 


(4.1.8) 


and the components of Eq. (4.1.2) become 

3v 3v 3v 

X ^ X ^ X 

— — + v + v • 

3t x 3x y 3y 



and 


3 V 

_J L 

3t 


3 V 

v + 

x 3x 


if = eg (t-t) + 



(4.1.9) 


Note that Eqs. (4.1.8) and (4.1.9) contain no terms for pressure (or modified 
p) which implies that (9p/9x). = ( 3 p/3 y ) = 0. This is only absolutely valid 
if, in fact, there is no convection. As the convection increases, Eqs. 
(4.1.8) and (4.1.4) become less valid; consequently, Eqs. (4.1.6) and (4.1.7) 
should be more appropriate to use over a wider range of convective conditions. 


4.2 Numerical Methods 

Previous investigators (26, 27-31) of natural convection in enclosures have 
all recast their approximations to Eq. (4.1.2) into a vorticity equation 
before seeking a numerical solution. Following the method of Wilkes and 
Churchill (26) the following procedure is employed: 

1. Eqs. (4.1.6) and (4.1.7) are differentiated with respect to y and x, 
respectively; 

2. the results are subtracted; and 

3. Eq. (4.1.1) is applied. 


This gives an equation in which pressure does not appear. 



Next, the vorticity is defined as 




3v 

_J t “ 
9x 



3y 


(4.2.1) 


(4.2.2) 
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and the stream function as 


v 


x 




Moreover, £ and 4 1 are related by (26, 32) 



(4.2.3) 


(4.2.4) 


By making use of Eq. (4.2.1), Eq. (4.2.2) becomes: 

If + v x ta + v y § ’ U + 

The governing equations are now summarized: 

Vorticity Equation - Eq. (4.2.5), 

Energy Equation - Eq. (4.1.3), 

Stream Function Equation - (Eq. (4.2.4), and 
Velocity Equation - Eq. (4.2.3). 

The relationships presented herein can be used to calculate the velocity and 
temperature fields in the bulk liquid provided the natural convection Is 
laminar. Since the system, comprising a S/L region encroaching upon a 
diminishing bulk liquid as depicted in Fig. 3.1* has not been studied for the 
limit of laminar stability, we examine related systems to predict the 
applicability of the governing equations. 


To determine whether flow is laminar or turbulent, the Rayleigh number is 
calculated. The form of the Rayleigh number and its critical value depend on 
the type of thermal boundary condition at x ** 0 (see Fig. 3.1) and the 
characteristic dimension selected for the system; several examples follow. 


Uniform Heat Flux at x = 0 

For a vertical surface of length Y, the Rayleigh number is defined by 


Ra = 


i| 

3 KS 

k T va 
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(4.2.6) 



in which q is the uniform heat flux applied at x = 0. 
convection is laminar (31)- 


If Ra < 


12 

10 , the 


In a rectangular enclosure of width X and length Y oriented with its length 
along the vertical, the Rayleigh number is defined as 

If Ra < 3 x 10^, the convection is laminar (31). 


Uniform Temperature at x ■ 0 

For a vertical surface of length Y, the Rayleigh number is 


Ra 


qBY 3 AT 

va 


(A. 2. 8) 


Here AT is the temperature difference between the vertical surface and the 

9 

liquid. If Ra <10, the convection is laminar (32). 

In a rectangular enclosure of width X and length Y oriented with its length 
along the vertical, the Rayleigh number is 


Ra 


q SX^AT 

va 


If Ra < 10 7 flow is laminar (33). 


(A. 2. 9) 


In order to estimate an appropriate value of q to be used in Eqs. (A. 2. 6) and 
(A. 2. 7), data on the cooling rate of Sn -15? Pb alloy before solidification 
were used. The alloy was contained in a mold of the dimensions Y «* 6.35 cm 
and X • 2.7 cm. Using the maximum cooling rate observed, 0.25 °C/S, the heat 
flux into the chill surface (i.e., q at x = 0) was computed by means of a 
simple energy applied to the alloy melt. 

A value of 25°C was used for AT in Eqs. (A. 2. 8) and (A. 2. 9). Since the 
temperature gradients within a convecting bulk liquid of a solidifying alloy 
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are small for the type of system considered, herein, the value of 25°C is 
quite conservative. 

The results of the calculations are summarized in Table 1. 


Table 1 

Summary of Calculations for Rayleigh Numbers 


Equation for 
Rayleigh 
Number 

Critical 

Rayleigh 

Number 

Calculated Rayleigh Numbers 

Sn-15S Pb 

A1-4S Cu 

Eq. (4.2.6) 

10 12 

10 6 

L 

5xl(T 

Eq. (4.2.7) 

3x10 6 

3x10** 

2x10 3 

Eq. (4.2.8) 

10 9 

10 6 

2xl0 5 

Eq. (4.2.9) 

10 7 

■H 

10 A 


In all cases, the Rayleigh numbers are well below those at which laminar 
convection would become unstable. Although the exact system is not duplicated, 
these numbers show that the assumption of laminar flow in the bulk liquid is 
valid. 

After solidification proceeds, the characteristic dimension (in the cases of 
Eqs. (4.2.7) and (4.2.9)) would decrease; also, q and/or AT would decrease. 
All of these factors reduce the Rayleigh numbers to values less than the 
calculated values shown in Table 1. 

Although Table 1 clearly shows that the convection in the bulk liquid is 
expected to be laminar, numerical calculations of previous workers indicate 
unstable results when the Grashof number exceeds a critical value. The 
unstable results have been attributed to the errors and approximations 
involved in the numerical schemes and should not be confused with a critical 
value of the Rayeligh number which indicates the laminar to turbulent 
transition. The Grashof number is defined as 
























Gr = 


(4.2.10) 


g bp Jut 

2 

w 

in which L is a characteristic dimension of the system (e.g., Y and X in Eqs. 
4. 2. 6-4. 2. 9) . The Rayleigh number and the Grashof number are related by 

Ra = Gr • Pr (4.2.11) 

in which Pr is the Prandtl number (v/a) of the bulk liquid. For the alloys 

considered in Table 1, the Prandtl numbers are approximately 0.02-0.03 so the 

2 4 

Grashof numbers would be in the range 5x10 to 2x10 for the Sn-15$ Pb alloy 
and 3x10^ to 7x10^ for the A1-4$ Cu alloy. 

Previous workers (8,26-30) who have sought numerical solutions to natural 

convection in enclosures report that difficulties in obtaining stable 

2 4 

solutions occur at Grashof numbers which exceed the range (10 - 10 ) to be 

encountered in this work. For example, in 1966 Wilkes and Churchill (26) 
report that they obtained numerical instabilities whenever the Grashof number 
exceeded 2x10"\ More recently (1981) Desai and Kim (28) reported stable 
numerical results for Grashof numbers in excess of 10^, and in a study of the 
bulk liquid convection above a solid-liquid region in a solidifying ingot, 
Ridder (8) achieved stable solutions for Grashof numbers up to approximately 
2x10 6 . 


4-7 



SECTION 5 
OPERATING GUIDE 


5.1 Running the Interactive Models 

In 1982 the models were converted to operate in interactive mode on a Digital 
Equipment Corporation VAX 11/780. The following discussion pertains to 
operating the model on a VAX, although only the starting procedure is system- 
dependent. The log-in procedure may vary from installation to installation 
and so is not described here. 

5.1.1 Starting the Solidification Model 

After logging in, the user initiates execution of the program by entering one 
command: MPS 3- The model identification page will appear on the screen 
shortly after this command is entered. 

5.1.2 Input Phase 

For each case, the program passes through three phases of operation. The 
first is the input phase during which the user defines the case to be run by 
specifying the case title, the alloy and the values of the model parameters. 
At each step the program clearly specifies the type of information it 
requires. The code is able to trap and recover from input errors, thus 
avoiding fatal FORTRAN I/O errors which would necessitate restarting the 
program. The input parameters are defined in a later section. Each time 
execution is initiated the built-in default case will appear in the model. If 
any parameter is changed during the input phase, the new value will remain in 
effect through subsequent cases in the same session or until it is changed 
again. Copies of the input phase displays are shown in the examples section. 

As described in reference 11, the input phase is controlled at a higher level 
so that the operator need only review the parameter classifications that are 
to be changed. After the case title is input the following selection appears 
on the screen: 
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ENTER P TO PROCEED TO CALCULATION 

1 TO CHANGE ALLOY 

2 TO REVIEW ALLOY DATA 

3 TO REVIEW OR CHANGE SOLIDIFICATION PROCESS PARAMETERS 
A TO REVIEW OR CHANGE PERMEABILITY MODEL PARAMETERS 

5 TO REVIEW OR CHANGE NUMERICAL METHODS CONTROL PARAMETERS 


If the operator selects one of 1 - 5 then the input phase proceeds directly to 
that classification. Upon exiting the classification, control returns to the 
higher level selection shown above. Thus the operator can move through the 

input at will. When he is satisfied with all input values, the "P" selection 

initiates the calculation. Some operational notes for each step of the input 
phase follow. 

o The case title is simply typed in; if none is desired, push 
RETURN. 

o For binary alloys, the solvent and solute are entered as 

(uppercase) chemical symbols with no leading or embedded 
blanks. The program then checks the alloy data base, and if 

the data for the requested alloy is found, it is displayed on 

the screen. If the alloy is not in the data base, an error 

message and a list of alloys in the data base are displayed, 
and the user is allowed to select another alloy. 

o The remaining parameters are grouped under several broad 

c lass if ic iat ions ; each group is a separate display page in the 
input phase. When a parameter is changed, the display is 

refreshed with the new value. The user can change any, all, or 
none of the parameters in a group before returning to the 

higher control level. 

5.1.3 Calculation Phase 

When the operator enters "P" during the input phase, the model begins the 

second stage of operation, the calculation phase. During this phase, the 

program determines the numerical solution to the equations describing the 

model subject to the conditions defined by the input parameters. 


After the input phase is complete, the model takes one step away from the 
chill and then prints the message shown below. 


5-2 



ENTER T TO DISPLAY TABULAR DATA, 

G TO DISPLAY GRAPHS, 

C TO CONTINUE THIS CASE, 

Q TO TERMINATE RUN, OR 
N TO PROCEED TO NEXT CASE. 

This choice of action is given to the user at each checkpoint. Entering a G 
or a T will cause the program to go to the output phase in which the user can 
select specific graphical or tabular output from the calculation up to that 
point. When the user leaves the output phase, the selection shown above is 
again presented. The user can switch freely between tabular and graphical 
output or he can continue the calculation, start a new case or terminate 
execution. If the continue (C) option is selected, the program will display 
the following: 

TIME SINCE BEGINNING OF SOL I Dl F I CAITON - 1.0000 SEC. 

ENTER THE TIME OF THE NEXT CHECK POINT. 

The checkpoint procedure puts the user in control of the progress of the 
calculation. The stopping places specified by the user as checkpoints allow 
him to inspect the progress of the calculation or display the current state of 
the S/L zone at any point in the solidification process. The stopping point 
is specified as a time simply by tying in the stopping time in seconds. For 
example, 

10 . 

The program then resumes the calculation and displays 

CONTINUING CALCULATION 

5.1. A Output Phase 

The three levels of operation during the output phase are illustrated in 
Figure 5.1. At the highest level the user can choose to display graphical 
output, to display tabular output, to terminate execution of the program, or 
to terminate the case by proceeding to the next case. 
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o At the graphical output level the user is given a choice of 
functions to plot and types of plot to use as shown in example 

5.4. Immediately after entry of the plot type, the plot will 

be generated on the screen*. When the cursor remains 

motionless at the upper left corner of the screen, the plot is 

complete. To leave the plot display and return to plot 

selection mode, enter a P. 

o At the tabular output level the display commences immediately 
after user selection of a function. Each time the screen in 
filled, the program waits for the user to enter a P before 
proceeding to the next display page or to the function 
selection page. On a typical 20 X 20 mesh a scalar function 
such as g. requires only one page while a vector function such 
as velocity requires three pages. 

o Note that it is possible to switch freely between graphical and 
tabular output at the highest control level and to repeat the 
plot of a given function with a variety of scaling choices. 

o After entering a Q to terminate execution or an N to proceed to 
the next case, the user can return to output mode only by re- 
running the case. 

5.2 Input Parameters 
5.2.1 Process Parameters 

The parameters in this group describe the conditions under which the casting 
is made including the mold geometry, the thermal conditions, and the strength 
and direction of the gravitational force. 


Parameter Description 

(Symbol ) 

Default 

Limits 

Mold half width 

X 

8 cm 

X > 0 

Mold height 

Y 

10 cm 

Y > 0 

Initial temperature of 
melt 

- - 

668 °C/sec 

greater 

than T lo 

Temperature of chill 

T ch i 1 1 

30 °C/sec 

T ch i 1 1 < T E 

Heat transfer coef- 
ficient at chill 

(h) 

0.1608 

h > 0 

Gravitational force 

(g) 

1 G 

g > 0 

Angle of g with x-axis 

' 

-90. 

-360 to 360 

* In some cases, the user 

is given the 

opportunity to override 

the 

computed scale with an 

input scale. 
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5.2.2 Permeability Model Parameters 

The permeability model is identical to that defined in Ref. 11. 


Parameters 

Symbol 

Default 

Limits 

Reference value of gamma 

Y o 

6 X 10-7 cm 2 Q 
0 

0 

A 

O 

>- 

Reference arm spacing 

d 

0 

1. micron 

d >0 
0 

Dendrite-arm spacing exponent 

q Y 

0 . 

O 

A 

of 

Solidification-time coefficient 

a d 

1. 

a .>0 
d 

Solidification-time exponent 

q d 

0. 

v°* 

y location of grain boundary 

- 

0. cm 

0. to I 

Y multiplier at grain boundary 

- 

1. 

0 


The default values were selected so that y = Y 0 with no grain boundary. 

5-2.3 Numerical Methods Control Parameters 

In the vast majority of cases the parameters in this group should not be 
changed from their default values, however operator control of these 


parameters could be useful 

in trouble-shooting 

• 


Parameter Description 

(Symbol) 

Default 

Limits 

Number of x-mesh points 

(n l ) 

10 

10 < n l < 40 

in the 1 iquid 




Number of x-mesh points 

<v 

10 

10 < N m < 40 

in the mush 




Number of x-mesh points 

(N $ ) 

10 

10 < N s < 40 

in the sol id 




Number of y-mesh points 

(N y ) 

10 

3 1 N y < 20 

Time step size 

(At) 

1 . 

At > 0. 


5.3 Alloy Data Base Structure 

The alloy data base contains alloy properties used by the solidification 
model. The use of a separate data base automatically accessed by the model 
frees the user from inputting relatively static alloy data with each change in 
mold contents, yet it allows complete flexibility in the range of alloys which 
can be processed. The phase diagram, density data and viscosity data for each 
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alloy are included as well as a reference block allowing notation of data 
sources. The data base is in card image format and can be modified by using 
any system text editor. 


CARD 

FORMAT 

CONTENTS 

1 

20A4 

Data base identifier 

2 

A i *,6X,A4,6X,2E10.3 

Solvent, solute, minimum C 
maximum C^ 

3 

20AA 

5 cards for notation of 
data sources 


ii 


5 

ii 


6 

ii 


7 

ii 


8 

4E10.3 

dC L /dT, k, C £ , T e 

9 

4E10.3 

d p L /dC L , p $ , p le , P SE 

10 

El 0 . 3 

y 

11 

7E10.3 

H L ,H S’ H E ,C PL’ C PSL ,C PSE ,C PE 

12 

4E10.3 

a L » 1 

13 

5E10.3 

a SL0 ,b SL0’ T HIL0 ,a SHI ,b SHI 

Cards 2 

through 13 are repeated for 

each alloy system. 
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5.** Sample Case (Cont.) 
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5. 4 Sample Case (Cont. 
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HHH UIM. AVERAGE COMPOSITION 
INTiROENDRITIC FLUID VELOCITY 




5.** Sample Case (Cont. 
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5.4 Sample Case (Cont. 
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6.1.2 Calculation Section 
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6.1.2 Calculation Section (Continued) 



CALCULATE X 


(k+1) 


CONVERGE? 


/ LESS ^ 
THAN MAXIMUM 
ITERATIONS 



EXIT WITH ERROR 
CONDITION 
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6.1.3 Input Section 














I I 
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6.2 Alphabetical Listing of Subroutines 

The list below contains only the modules that were written specifically for 
the solidification model. VAX/VUS routines and Tektronix routines called by 
the model are listed in Ref. 11. 


Ca 1 cu 1 at i on Modules 


NAME 

CNTRL 


CROSS 
DNSTY 
ENGY1 ,2,3 

FINSLD 

GLlfclT 

HOFT 

HTEVAL 

IN1T 

LQD 

MASS2 

MSH 

PENTA 

PERM 

PMLQD 

PMMSH 

PMSLD 

SLD 


FUNCTION 

Tests for and controls change of calculation phase, for example 
from all-liquid conditions to the beginning of the initial 
transient. 

Locates the position where temperature has an input value. 
Calculates functions associated with density. 

Evaluate Q and the Jacobian elements for the energy equation in 
the liquid, mush and solid, respectively. 

Calculates macrosegregation in the final solid. 

Interpolates volume fraction liquid. 

Evaluates enthalpy as a function of temperature. 

Evaluates Q and the Jacobian elements for the enthalpy- 

temperature relationship in each of the three regions. 
Initialization routine. 

Controls the calculation of information for the bulk liquid 
region during the Newton-Raphson iterations. 

Evaluates Q and the Jacobian elements for the total-mass 

conservation equation in the mushy zone. 

Controls the calculation of information for the mushy zone during 
the Newton-Raphson iterations. 

Solves penta-diagonal matrix equations. 

Evaluates the permeability function. 

Performs partitioned matrix calculations for elements from the 

bulk 1 iquid region. 

Performs partitioned matrix calculations for elements from the 

mush zone. 

Performs partitioned matrix calculations for elements from the 

solid region. 

Controls the calculation of information for the solid region 

during the Newton-Raphson iterations. 
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SOLU 


Evaluates Q and the Jacobian elements for the solute-mass 
conservation equation in the mushy zone. 

SPLIT Called by CNTRL to separate the temperature array into arrays 

over two regions. 

STEP Controls the Newton-Raphson iterations. 

STORE Stores arrays from the previous time prior to taking another time 

step. 

STRECH Interpolates a static array over a growing region. 

THCOND Evaluates the thermal conductivity. 

TRI2,3 Solve tridiagonal matrix equations. 

VLCTY Evaluates the velocity from Darcy's Law. 

XYCALC Calculates some of the coefficients of the coordinate 

transformation. 


Input Modules 

All input modules are the same as those listed in Ref. 11. 


Graphics Modules 


NAME 

AXES 

CMPVEC 

CTRBOX 

CTRDRW 

DRWVEC 

FRACEL 

GPBLK 

GPHCON 

GPHG 

HPROFS 

INGOT 

ISOTH 

PRFINT 

SCALE 

SETPLT 

VPROFS 

XL ABEL 


FUNCTION 

Draws and labels axes. 

Plots composition contours and velocity vectors. 

Controls contouring within a grid cell. 

Draws contour lines. 

Draws arrows for vector plots. 

Plots contours of fraction eutectic and fraction liquid. 
Puts parameter block on plots. 

Controls graphical output. 

Puts gravity vector on plots. 

Draws horizontal profiles. 

Sets up axes and labels for plots on whole ingot. 

Plots isotherms throughout casting. 

Interpolates profiles between grid points. 

Calculates neat scaling. 

Initial izes plot. 

Draws vertical profiles. 

Sets up label for horizontal axis. 
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6.3 Key Program 

Symbol s 

✓ 

PROGRAM 

COMMON 


SYMBOL 

BLOCK 

DEFINITION 

ACS 

/AUX/ 

‘S 

CL 

/AUX/ 

C L 

DX 11 , DX 1 2 , DX 1 3 

/MESH/ 


DYE 

/MESH/ 

An 

GL 

/STATE/. ' 

9l 

GE 

/AUX/ 

£e 

H1,H2,H3 

/STATE/ 

pH in the liquid, 

N 11 , N 1 2 , N 1 3 

/MESH/ 

n l ,n m ,n s 

NJ 

/MESH/ 

N y 

P2 

/STATE/ 

p 

RL2 

/AUX/ 

p L 

T1,T2,T3 

/STATE/ 

T in the liquid, 

V 

/AUX/ 

V 

XC 

/MESH/ 

X 

X 1 1 , X 12 , X 13 

/MESH/ 

? L ,? M’ ? S 

XL,XE 

/MESH/ 

and at each 

YE 

/MESH/ 

n 

YY 

/CTSFM/ 

Y 


solid 


solid 
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APPENDIX A 
THERMAL PROPERTIES 


A. 1 Enthalpy - Temperature Coefficients 

The following values were derived from Appendix B which is an analysis of the 
enthalpies of the solidifying alloys. 



A1 -k . 5 wt. pet. Cu 

Sn-15 wt. pet. Pb 

PL 

1.059 (J/gm°C) 

0.237 ( J/gm°C) 

PSL 

1.310 ( J/gm°C) 

0.250 ( J/gm°C) 

PSE 

1.033 ( J/gm°C) 

0.2^3 ( J/gm°C) 

PE 

0.837 ( J/gm°C) 

0.161 ( J/gm°C) 

: 

391.6 (J/gm) 

54.27 (J/gm) 

s* 

0. (J/gm) 

0. (J/gm) 

E 

-277.8 (J/gm) 

-13.39 (J/gm) 


A. 2 Thermal Conductivity of Al-4.5 wt. pet. 
(All thermal conductivities are in units J/cm 

solid: k T$ - -2.09x10 -3 (T-T E ) + 0.709 

k T$ = 5.27x10"**' (T-T e ) + 1.286 

liquid: k JL = A. 18x10"** (T-T E ) + 0.867 

S/L: k T = g $ k TS + g L k TL 


Cu (Ref. 34-36) 

. s . °C . ) 

if 27°C < T < 327°C 
if 327°C < T < 6*»5°C 
if T > 5W°C 


A. 3 Thermal Conductivity of Sn-15 wt. pet. Cu (Ref. 37-39) 


(All thermal 

conductivities are in units J/cm.s.°C.) 

sol id: 

k T ’ % k TS + 9 E k TE 

with 

k JS = 39x10"** (T-T E ) + 0.599 if T < 183°C 


k TE => -2.30x10"** (T-T £ ) + 0.^26 if T < l83°C 

1 iquid: 

k JL - 2.07x10"** (T-T e ) + 0.292 if T > l83°C 

S/L: 

k T * 9 S k TS * Vti 

with 

k JS = -2.70x10"** (T-T e ) + 0.609 if T > l83°C 
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APPENDIX B 

SOLIDIFICATION ENTHALPIES AND THE EFFECT OF 
MACROSEGREGATION ON ENTHALPY 


B. 1 The Enthalpy of the Primary Solid 

The solid is assumed to be a regular solution so that its enthalpy (for a 
binary system) is given by 


X,H° S + X 2 H° s + y,X 2 
X 1 M 1 + X 2 M 2 


(B. 1 ) 


where 

H,. = enthalpy of the solid (cal/g) 

X^,X^ = mole fractions of components 1 and 2, respectively, 

H ?S’ H 2S = mo ^ ar enthalpies of the pure components as solids (cal/mole), 

« s ■ constant which depends upon the heat of mixing for the solid 

(cal/mole), and 

Mj,M 2 ■ atomic weights (g/mole). 


The enthalpies of the pure components, themselves, are temperature dependent 
and are given by 


H? - a, (T-298) + -j (T 2 -298 2 ) + c. (298" 1 -T" 1 ) 


(B.2) 


where 

T = absolute temperature (K), 

and aj, b. and c^ are empirical constants for the specific heats. In 
particular, the specific heats are expressed as 

C p . = a. + b.T + c.T" 2 (B.3) 


with C . in (cal/mole-K). 
P' 
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B.2 The Enthalpy of the Liquid 

The liquid is assumed to be a regular solution so that its enthalpy (for . a 
binary) is given by 


x, h “l * x 2 h£ l ♦ « l X,X 2 
x 1 'm 1 + x~M 2 


(B.4) 


* enthalpy of the liquid (cal/g), 

= mo ^ ar enthalpies of the pure components as liquids (cal/mole), 


«= constant which depends upon the heat of mixing for the liquid 
(cal/mole). 


The enthalpies of the pure components are 


h ° l = (H°i m + AH fj) + / < d j + e s T + 9 j T ” 2 ) dT 

mi 


(B.5) 


where 

H? m ■ enthalpy of the pure solid 1 at its melting point (cal/mole) 

AH^.. ** heat of solidification for pure i (cal/mole), 

T . = melting point of pure i (K), and 

d.,e.,g. » empirical constants for the molar specific heat of liquid i. 

I 

The integration of Eq. (B.5) yields: 


H ?L " <H ?n, * iH fi> + (T2 - T ml> 


+ g,(T"! - T’ 1 ) 
i mi 


(B.6) 
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B.3 The Enthalpy of the Eutectic Constituent 
Two types of eutectic are considered: 

1. the eutectic is a mixture of a terminal and an intermetal 1 ic 
phase (e.g., in Al-Cu alloys), and 

2. the eutectic is a mixture of two terminal phases (e.g., in Sn-Pb 
al loys). 

When the eutectic comprises a terminal phase and an intermetal 1 ic phase, the 
enthalpy is 

H E " f | H | + (1_f l )H a (B.7) 

in which 

= enthalpy of the eutectic mixture (cal/g), 

Hj = enthalpy of the intermetal 1 ic phase (cal/g), 

= enthalpy of the terminal phase (cal/g), and 
fj = weight fraction of the intermetal 1 ic phase in the eutectic mixture. 

The enthalpy of the intermeta 1 1 ic phase is estimated by using the Neumann-Kopp 
rule (AO) for its molar specific heat; accordingly 


C p| - mC pi + nC p2 


(B.8) 


and 



AH 29 g + mH“ + nH 

mMj + nM 2 


0 

2 


(B.9) 


Here, 

C p | =» molar heat capacity of the intermetal 1 ic phase (cal/mole-K) , 
m,n = stoichiometric proportions in the intermetal 1 ic phase, and 
AH 298 “ heat ^ orma ^ion at 298K. 


in Eq. (B.7), H q is evaluated by Eq. (B.l) in which = H a and X 2 is the mole 
fraction of solute in the terminal phase within the eutectic mixture. 
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When the eutectic comprises two terminal phases, the enthalpy is 


f a H a 


+ 


(B. 10) 


where 

H a ,Hg = enthalpies of the terminal solid solutions making up the eutectic 
(cal/g), and 

f Q ,fg = weight fractions of the terminal solid solutions in the eutectic. 


In Eq. (B.10), H q and are evaluated by Eq. (B.l) with X 2 as the mole 
fraction of solute in the respective terminal phases in the eutectic mixture. 


B.A The Enthalpy of Sol idifyinq Al loys ' 

The enthalpy of the alloy when it comprises both solid and liquid is 

PH = ( 1 - 9 L )P $ H S + g L P L H L T e < T < T lq (B.l 1 ) 

where, 

Pg,P L = mass densities of the solid and liquid phases (g/cm^) 
respectively, 

g L - volume fraction of liquid, and 

P = ( 1 ~gj_ ) P 5 + 9l p l> avera 9 e density of the solid-liquid mixture. 

Upon complete solidification, the alloy comprises the eutectic constituent and 
the primary solid. Thus 

J 

PH - (l-g E )P s H $ + 9 E P SE H E T < T e (B.12) 

where, 

g F = volume fraction of eutectic, and 

^ 3 

P^ E = mass density of the eutectic (g/cm ). 

In this case, the average density of the solidified alloy Is 

p “ o-g E )p s + 9 e p se • 
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B.5 The Enthalpy of A1-^.5% Cu and Sn-15% Pb Alloys During Solidification 
During solidification and after solidification as the alloy cools, the 
enthalpy is a function of the local macrosegregation since enthalpy depends 
upon the relative amounts of the constituents and the composition of the 
constituents. Specifically, in the S/L temperature range the relative amounts 
of the phases as well as the local average composition of the dendritic solid 
depends upon the specific solidification conditions at that locale, and after 
solidification the local fraction of the eutectic constituent (and hence local 
average composition) can vary substantially from point-to-point within the 
costing or ingot. Hence, enthalpy can not be treated as a function of 

temperature only. Preliminary calculations showed, indeed, that enthalpy 
depends strongly on the degree of macrosegregation at a specific location. 
However, the effect of the average composition of the primary dendritic solid 

on the enthalpy of the solid can be neglected. This simplifies the computation 

of enthalpy in that the enthalpy of the alloy comprises the enthalpies of the 
liquid, solid and eutectic constituents in accordance with their relative 

amounts regardless of the average composition of the primary solid. 

Figure B1 shows the results of computing enthalpy for the three constituents 
in Al-i».5$ Cu alloy. For the liquid-phase there is a unique curve since the 
composition of the liquid above and below the liquidus temperature is a 
function only of temperature regardless of the local flow and macrosegregation 
conditions. For the primary solid two curves are shown, each corresponding to 
a very different macrosegregation condition. in one case, the net flow of 
interdendr i t ic liquid has been strongly positive resulting a volume fraction 
of eutectic of 0.20 and an average composition after solidification of 9.5$ 
Cu. In the other case, the volume fraction of eutectic is 0.066 and the 
average composition is only 3.^5$ Cu. Thus, one curve is for strong positive 
segregation, and the other is for negative segregation yet the enthalpies are 
almost equal. Therefore, the enthalpy of the primary solid is assumed to be 
only a function of temperature. Finally, there is only one curve for the 
eutectic constituent since it is assumed that the composition of the eutectic 
which forms at the end of solidification is unique. 

The enthalpy of the alloy in the S/L temperature range (T £ ± T ± T L ) is given 
by the weighted average of the enthalpies of the liquid and primary solid 
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Figure B1 - Enthalpy of the liquid and solid phases and eutectic 
constitutent for Al-4.5% Cu 
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ENTHALPY, cal 



Figure B2 - Enthalpy of liquid and solid phases and eutectic 
constituent for Sn - 15% Pb 

Solid line for th£ primary solid corresponds to positive 
segregation wi th C = 21.5.% Pb and the broken line to 
negative segregation with C g = 13-2% Pb. 
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constituents, with the relative amounts of the constituents determined by the 
local macrosegregation conditions. Likewise for T _< T^, the enthalpy of the 
alloy depends upon the relative amounts of the primary solid constituent and 
the eutectic constituent. Thus, the calculation of enthalpy is a dynamic 
calculation which is done as part of the entire macrosegregation calculation. 

Figure B.2 shows the enthalpy curves for Sn-15$ Pb alloy. Again the two 
curves for the primary solid are very close to each other so the enthalpy of 
the primary solid is approximated as a unique function of temperature. 
Calculations for Figures B.1 and B.2 were performed using the thermodynamic 
properties listed in Table B.1. Linear regression analyses on many 
coordinates along each of the curves shown in Figures B.1 and B.2 were 
performed to obtain the values listed in Appendix A.1. 



Table B.l 


Thermodynamic properties of AI-4.5% Cu 
and Sn-151 Pb Alloys 


Alloy 

A1-4.5* Cu 

Sn-15% Pb 

Refer- 

ence 

Components 

A1 

Cu 

Sn 

Pb 

a. (cal/mol-K) 

4.94 

5.41 

4.42 

5.63 

41 


mmm 

0.00150 

0.0063 

■BBSS 

41 

c. (cal-K/mol) 

0 

0 

0 

0 

41 


7.0 

7*5 

8.29 

7.75 

41 

f. (cal/mol-K^) 

0 

0 

-0.0022 

-0.00074 

41 

g. (cal-K/mol) 

0 

0 

0 

0 

41 

T j (K) 
m 1 

933 

1,356 

505.1 

600.6 

42,44 

AH^.. (cal /mol) 

2,550 

3,120 

1,680 

1,147 

42,44 

a $ (cal /mol) 

-19,100 

3,244 (1) 

43 

« L (cal /mol) 

-5,550 

1,319 

43 

m 

2.13 

— 

43 

n 

1 

— 

43 

AH 2 gg (cal /mol) 

-9,600 

— 

43 

f l 

0.584 


1 

1 

X a 

0.0248 

0.985 

43 

X B 

— 

0.71 

43 


(1) This is approximate since data for Pb-rich alloys are used. 
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